##########
#Replication Code
#Project: When Green Trade Backfires: The Uneven Effects of the EU Timber Regulation on Deforestation 
#Author: Stefano Jud and Quynh Nguyen
#Date: January 19, 2026
##########

########Preamble----

rm(list = ls())



if (!require("pacman")) install.packages("pacman")

pacman::p_load(devtools, tidyverse,fect,gsynth,fixest,rnaturalearth,ggfixest,ggpubr,
               marginaleffects,stargazer,plotrix,lfe,did,zoo)

devtools::install_github("synth-inference/synthdid")

library(synthdid)

#setwd() #Enter string with your working directory

##### Load Main Data ----

df_no_EU <- read.csv('data_JudNguyen_EUTR.csv')


##################################
#####FIGURES MAIN MANUSCRIPT######
##################################

#####Figure 1----- 

treat_countries <- df_no_EU %>% filter(pre_EUTR_eu_top_imp_weight == 1 & !is.na(deforest_30)) %>% distinct(iso3c) %>% as.matrix()
control_countries <- df_no_EU %>% filter(pre_EUTR_eu_top_imp_weight == 0 & !is.na(deforest_30)) %>% distinct(iso3c) %>% as.matrix()

world_coordinates <- ne_countries(scale = "medium", returnclass = "sf") %>%
  mutate(treated = if_else(iso_a3 %in% treat_countries,'Treated Countries',
                           if_else(iso_a3 %in% control_countries,'Control Countries','Not in Sample'))) %>%
  filter(su_a3 != 'ATA') 

ggplot() + 
  geom_sf(data=world_coordinates,aes(fill = treated)) + 
  scale_fill_manual(breaks = c("Treated Countries", "Control Countries", "Not in Sample"), 
                    values=c( "#FFA500","#006696", "grey50")) +
  labs(fill = 'Treatment Group') +
  theme_void() 

ggsave('Output//Fig1.jpeg',width=16,height=6,dpi = 400)

#####Figure 2----- 
mean_deforest <- df_no_EU %>% drop_na(pre_EUTR_eu_top_imp_weight) %>% 
  group_by(pre_EUTR_eu_top_imp_weight,year) %>% 
  summarise(mean_deforest = mean(log_deforest_30,na.rm = T),
            mean_deforest_rate = mean(deforest_percent_30,na.rm = T)) %>% ungroup()


ggplot(mean_deforest %>% filter(year < 2020),
       aes(x=year,y=mean_deforest,shape = as.factor(pre_EUTR_eu_top_imp_weight),linetype = as.factor(pre_EUTR_eu_top_imp_weight))) + 
  geom_point(size = 3) + geom_line(linewidth = 0.8) +  
  scale_shape(labels = c("Control Countries","Treated Countries")) +
  scale_linetype(labels = c("Control Countries","Treated Countries")) +
  labs(linetype = 'Treatment Group',shape = 'Treatment Group') + 
  scale_x_continuous(breaks=seq(2001,2019,1)) +
  geom_vline(xintercept = 2010,linetype = 'dashed') + 
  xlab('Year') + ylab('Mean Log Deforestation (in ha)') +
  theme_bw()

ggsave('Output//Fig2.jpeg',width=12,height=6,dpi = 400)

#####Figure 3----- 
ALL_es <- feols(log_deforest_30  ~ i(year, pre_EUTR_eu_top_imp_weight, 2009) |                    ## Other controls
                  year + iso3c,                             ## FEs
                cluster = ~iso3c,                          ## Clustered SEs
                data = df_no_EU %>% filter(year < 2020))

jpeg(
  filename = "Output//Fig3.jpeg",
  width = 10,height = 7,units = "in",res = 400)

iplot(ALL_es,main = '',xlab = 'Year',ylab = 'ATT of EUTR on Log Deforestation',
      pt.join = F,lwd = 1.5,cex =2)
dev.off()

aggr_es(ALL_es, aggregation = "mean") #Describes aggregate treatment effect

#####Figure 4----- 

eu_share_change <- feols(eu_share_weight  ~ i(year, pre_EUTR_eu_top_imp_weight, 2009)|                    ## Other controls
                           year + iso3c,                             ## FEs
                         cluster = ~iso3c,                          ## Clustered SEs
                         data = df_no_EU %>% filter(year < 2020))

jpeg(
  filename = "Output//Fig4.jpeg",
  width = 10,height = 7,units = "in",res = 400)
iplot(eu_share_change,main = '',xlab = 'Year',ylab = 'Difference Timber Export Share to EU',
      pt.join = F,lwd = 1.5,cex =2)
dev.off()

aggr_es(eu_share_change,period = 'post',aggregation = 'mean')

#####Figure 5----- 
###Subset data to before Year < 2020
df_no_EU_2019 <- df_no_EU %>% filter(year < 2020)

###Model with Political Corruption (Panel A)
#Run Model
eu_share_change_het1 <- feols(eu_share_weight ~ i(year, pre_EUTR_eu_top_imp_weight, 2009)*pre_EUTR_vdem_pol_corrupt  |                    ## Other controls
                                year + iso3c,                             ## FEs
                              cluster = ~iso3c,                          ## Clustered SEs
                              data = df_no_EU_2019)

#Make Figure
corruption_colors <- c("#2ca02c", "#ff7f0e", "#d62728")  # Green → Orange → Red

het1_plot <- plot_slopes(
  eu_share_change_het1,
  variables = 'pre_EUTR_eu_top_imp_weight',
  condition = list(
    year = unique(df_no_EU_2019$year),
    pre_EUTR_vdem_pol_corrupt = c(0.86, 0.65, 0.2)  # Now: High → Medium → Low
  )
) +
  geom_hline(yintercept = 0, linetype = 'dashed') +
  labs(
    title = 'A) Political Corruption',
    x = 'Year',
    y = 'Difference in Share of Timber Exports to EU'
  ) +
  theme_bw() +
  ylim(c(-50, 35)) +
  theme(
    legend.position = "none",
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank()
  ) +
  scale_color_manual(values = corruption_colors) +
  scale_x_discrete(breaks = seq(min(df_no_EU_2019$year), max(df_no_EU_2019$year), by = 3)) +
  facet_wrap(
    ~ factor(pre_EUTR_vdem_pol_corrupt, levels = c(0.86, 0.65, 0.2)),  # Force order inside plotting
    ncol = 3,
    labeller = as_labeller(c(
      `0.86` = "High corruption score (0.86)",
      `0.65` = "Median corruption score (0.65)",
      `0.2`  = "Low corruption score (0.2)"
    ))
  )

#Estimate pooled effect
avg_slopes(eu_share_change_het1, 
                          variables = "pre_EUTR_eu_top_imp_weight", 
                          newdata = datagrid(
                            year = seq(2010,2019,1),
                            pre_EUTR_vdem_pol_corrupt = c(0.2, 0.65, 0.86)
                          ),by = "pre_EUTR_vdem_pol_corrupt")

###Model with Rule of Law (Panel B)
#Run Model
eu_share_change_het2 <- feols(eu_share_weight ~ i(year, pre_EUTR_eu_top_imp_weight, 2009)*pre_EUTR_vdem_ruleoflaw  |                    ## Other controls
                                year + iso3c,                             ## FEs
                              cluster = ~iso3c,                          ## Clustered SEs
                              data = df_no_EU_2019)

#Make Figure
rol_colors <- c("#d62728","#ff7f0e", "#2ca02c")  # Green → Orange → Red

het2_plot <- plot_slopes(eu_share_change_het2,variables = 'pre_EUTR_eu_top_imp_weight' , condition = list(
  year = unique(df_no_EU_2019$year),  # All unique years from the dataset
  pre_EUTR_vdem_ruleoflaw = c(0.13, 0.44, 0.83)  # Values of interest
)) + geom_hline(yintercept = 0, linetype = 'dashed') +
  labs(title = 'B) Rule of Law' ,x = 'Year', y = 'Difference in Share of Timber Exports to EU') +
  theme_bw()  +ylim(c(-50,35)) +
  theme(legend.position = "none") +  # Move legend for better readability
  scale_color_manual(values = rol_colors) +  # Adjust colors for clarity
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) + 
  scale_x_discrete(breaks = seq(min(df_no_EU_2019$year), max(df_no_EU_2019$year), by = 3)) +  # Show only every third year
  facet_wrap(~pre_EUTR_vdem_ruleoflaw, ncol = 3, labeller = as_labeller(c(
    `0.13` = "Low Rule of Law score (0.13)",
    `0.44` = "Median Rule of Law score (0.44)",
    `0.83` = "High Rule of Law score (0.83)"
  )))  # Separates estimates into clear facets

#Estimate pooled effect
avg_slopes(eu_share_change_het2, 
                          variables = "pre_EUTR_eu_top_imp_weight", 
                          newdata = datagrid(
                            year = seq(2010,2019,1),
                            pre_EUTR_vdem_ruleoflaw = c(0.13, 0.44, 0.83)
                          ),by = "pre_EUTR_vdem_ruleoflaw")

###Model with State Capacity (Panel C)
#Run Model
eu_share_change_het3 <- feols(eu_share_weight ~ i(year, pre_EUTR_eu_top_imp_weight, 2009)*pre_EUTR_state_cap  |                    ## Other controls
                                year + iso3c,                             ## FEs
                              cluster = ~iso3c,                          ## Clustered SEs
                              data = df_no_EU_2019)

#Make Figure
rol_colors <- c("#d62728","#ff7f0e", "#2ca02c") 

het3_plot <- plot_slopes(eu_share_change_het3,variables = 'pre_EUTR_eu_top_imp_weight' , condition = list(
  year = unique(df_no_EU_2019$year),  # All unique years from the dataset
  pre_EUTR_state_cap = c(-0.44, 0.23, 0.9)  # Values of interest
)) + geom_hline(yintercept = 0, linetype = 'dashed') +
  labs(title = 'C) State Capacity' ,x = 'Year', y = 'Difference in Share of Timber Exports to EU') +
  theme_bw()  + ylim(c(-50,35)) + 
  theme(legend.position = "none") +  # Move legend for better readability
  scale_color_manual(values = rol_colors) +  # Adjust colors for clarity
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) + 
  scale_x_discrete(breaks = seq(min(df_no_EU_2019$year), max(df_no_EU_2019$year), by = 3)) +  # Show only every third year
  facet_wrap(~pre_EUTR_state_cap, ncol = 3, labeller = as_labeller(c(
    `-0.44` = "Low State Capacity score (-0.44)",
    `0.23` = "Median State Capacity score (0.23)",
    `0.9` = "High State Capacity score (0.9)"
  )))  # Separates estimates into clear facets

#Estimate pooled effect
avg_slopes(eu_share_change_het3, 
                          variables = "pre_EUTR_eu_top_imp_weight", 
                          newdata = datagrid(
                            year = seq(2010,2019,1),
                            pre_EUTR_state_cap = c(-0.44, 0.23, 0.9)
                          ),by = "pre_EUTR_state_cap")

###Add all three Figures together
ggpubr::ggarrange(het1_plot,het2_plot,het3_plot,nrow = 3,ncol = 1)

ggsave('Output//Fig5.jpeg',width = 18,height = 18,dpi = 400)


#####Figure 6----- 

###Model with Rule of Law (Panel A)
#Run Model
china_share_change_het1 <- feols(china_share_weight ~ i(year, pre_EUTR_eu_top_imp_weight, 2009)*pre_EUTR_vdem_pol_corrupt |                    ## Other controls
                                   year + iso3c,                             ## FEs
                                 cluster = ~iso3c,                          ## Clustered SEs
                                 data = df_no_EU_2019)

#Make Figure
het1_plot_chn <- plot_slopes(
  china_share_change_het1,
  variables = 'pre_EUTR_eu_top_imp_weight',
  condition = list(
    year = unique(df_no_EU_2019$year),
    pre_EUTR_vdem_pol_corrupt = c(0.86, 0.65, 0.2)  # Now: High → Medium → Low
  )
) +
  geom_hline(yintercept = 0, linetype = 'dashed') +
  labs(
    title = 'A) Political Corruption',
    x = 'Year',
    y = 'Difference in Share of Timber Exports to China'
  ) +
  theme_bw() +
  ylim(c(-25,40)) +
  theme(
    legend.position = "none",
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank()
  ) +
  scale_color_manual(values = corruption_colors) +
  scale_x_discrete(breaks = seq(min(df_no_EU_2019$year), max(df_no_EU_2019$year), by = 3)) +
  facet_wrap(
    ~ factor(pre_EUTR_vdem_pol_corrupt, levels = c(0.86, 0.65, 0.2)),  # Force order inside plotting
    ncol = 3,
    labeller = as_labeller(c(
      `0.86` = "High corruption score (0.86)",
      `0.65` = "Median corruption score (0.65)",
      `0.2`  = "Low corruption score (0.2)"
    ))
  )

#Estimate pooled effect
avg_slopes(china_share_change_het1, 
                           variables = "pre_EUTR_eu_top_imp_weight", 
                           newdata = datagrid(
                             year = seq(2010,2019,1),
                             pre_EUTR_vdem_pol_corrupt = c(0.2, 0.65, 0.86)
                           ),by = "pre_EUTR_vdem_pol_corrupt")

###Model with Rule of Law (Panel B)
#Run Model
china_share_change_het2 <- feols(china_share_weight ~ i(year, pre_EUTR_eu_top_imp_weight, 2009)*pre_EUTR_vdem_ruleoflaw  |                    ## Other controls
                                   year + iso3c,                             ## FEs
                                 cluster = ~iso3c,                          ## Clustered SEs
                                 data = df_no_EU_2019)

#Make Figure
het2_plot_chn <- plot_slopes(china_share_change_het2,variables = 'pre_EUTR_eu_top_imp_weight' , condition = list(
  year = unique(df_no_EU_2019$year),  # All unique years from the dataset
  pre_EUTR_vdem_ruleoflaw = c(0.13, 0.44, 0.83)  # Values of interest
)) + geom_hline(yintercept = 0, linetype = 'dashed') +
  labs(title = 'B) Rule of Law' ,x = 'Year', y = 'Difference in Share of Timber Exports to China') +
  theme_bw()  +ylim(c(-25,40)) +
  theme(legend.position = "none") +  # Move legend for better readability
  scale_color_manual(values = rol_colors) +  # Adjust colors for clarity
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) + 
  scale_x_discrete(breaks = seq(min(df_no_EU_2019$year), max(df_no_EU_2019$year), by = 3)) +  # Show only every third year
  facet_wrap(~pre_EUTR_vdem_ruleoflaw, ncol = 3, labeller = as_labeller(c(
    `0.13` = "Low Rule of Law score (0.13)",
    `0.44` = "Median Rule of Law score (0.44)",
    `0.83` = "High Rule of Law score (0.83)"
  )))  # Separates estimates into clear facets

#Estimate pooled effect
avg_slopes(china_share_change_het2, 
                             variables = "pre_EUTR_eu_top_imp_weight", 
                             newdata = datagrid(
                               year = seq(2010,2019,1),
                               pre_EUTR_vdem_ruleoflaw = c(0.13, 0.44, 0.83)
                             ),by = "pre_EUTR_vdem_ruleoflaw")

###Model with State Capacity (Panel C)
#Run Model
china_share_change_het3 <- feols(china_share_weight ~ i(year, pre_EUTR_eu_top_imp_weight, 2009)*pre_EUTR_state_cap  |                    ## Other controls
                                   year + iso3c,                             ## FEs
                                 cluster = ~iso3c,                          ## Clustered SEs
                                 data = df_no_EU_2019)

#Make Figure
het3_plot_china <- plot_slopes(china_share_change_het3,variables = 'pre_EUTR_eu_top_imp_weight' , condition = list(
  year = unique(df_no_EU_2019$year),  # All unique years from the dataset
  pre_EUTR_state_cap = c(-0.44, 0.23, 0.9)  # Values of interest
)) + geom_hline(yintercept = 0, linetype = 'dashed') +
  labs(title = 'C) State Capacity' ,x = 'Year', y = 'Difference in Share of Timber Exports to China') +
  theme_bw()  + ylim(c(-25,40)) +
  theme(legend.position = "none") +  # Move legend for better readability
  scale_color_manual(values = rol_colors) +  # Adjust colors for clarity
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) + 
  scale_x_discrete(breaks = seq(min(df_no_EU_2019$year), max(df_no_EU_2019$year), by = 3)) +  # Show only every third year
  facet_wrap(~pre_EUTR_state_cap, ncol = 3, labeller = as_labeller(c(
    `-0.44` = "Low State Capacity score (-0.44)",
    `0.23` = "Median State Capacity score (0.23)",
    `0.9` = "High State Capacity score (0.9)"
  )))  # Separates estimates into clear facets

#Estimate pooled effect
avg_slopes(china_share_change_het3, 
                             variables = "pre_EUTR_eu_top_imp_weight", 
                             newdata = datagrid(
                               year = seq(2010,2019,1),
                               pre_EUTR_state_cap = c(-0.44, 0.23, 0.9)
                             ),by = "pre_EUTR_state_cap")

###Add all three Figures together
ggpubr::ggarrange(het1_plot_chn,het2_plot_chn,het3_plot_china,nrow = 3,ncol = 1)

ggsave('Output//Fig6.jpeg',width = 18,height = 18,dpi = 400)

#####Figure 7----- 

###Panel A: Commodity-driven deforestation
commodity <- feols(log(commodity_driven+1)  ~ i(year, pre_EUTR_eu_top_imp_weight, 2009)  |                    ## Other controls
                     year + iso3c,                             ## FEs
                   cluster = ~iso3c,                          ## Clustered SEs
                   data = df_no_EU_2019)
plot_commodity <- ggiplot(commodity,geom_style = 'ribbon') + ggtitle('A) Commodity-Driven Deforestation') + 
  xlab('Year') + ylab('Difference in Commodity-Driven Deforestation')

aggr_es(commodity) #Pooled effect

###Panel B: Forestry-driven deforestation
forestry <- feols(log(forestry_driven+1)  ~ i(year, pre_EUTR_eu_top_imp_weight, 2009)  |                    ## Other controls
                    year + iso3c,                             ## FEs
                  cluster = ~iso3c,                          ## Clustered SEs
                  data = df_no_EU_2019  )

plot_forestry <- ggiplot(forestry,geom_style = 'ribbon') + ggtitle('B) Forestry-Driven Deforestation') + 
  xlab('Year') + ylab('Difference in Forestry-Driven Deforestation')

aggr_es(forestry) #Pooled effect

###Panel C: Shifting Agriculture-driven deforestation
shifting_agr <- feols(log(shifting_agr_driven+1)  ~ i(year, pre_EUTR_eu_top_imp_weight, 2009)  |                    ## Other controls
                        year + iso3c,                             ## FEs
                      cluster = ~iso3c,                          ## Clustered SEs
                      data = df_no_EU_2019 )

plot_shifting_agr <- ggiplot(shifting_agr,geom_style = 'ribbon') + ggtitle('C) Shifting Agriculture-Driven Deforestation') + 
  xlab('Year') + ylab('Difference in Shifting Agriculture-Driven Deforestation')

aggr_es(shifting_agr) #Pooled effect

###Panel D: Urbanization-driven deforestation
urban <- feols(log(urban_driven+1)  ~ i(year, pre_EUTR_eu_top_imp_weight, 2009)  |                    ## Other controls
                 year + iso3c,                             ## FEs
               cluster = ~iso3c,                          ## Clustered SEs
               data = df_no_EU_2019 )


plot_urban <- ggiplot(urban,geom_style = 'ribbon') + ggtitle('D) Urbanization-Driven Deforestation') + 
  xlab('Year') + ylab('Difference in Urbanization-Driven Deforestation')

aggr_es(urban) #Pooled effect

###Compile Figures
ggpubr::ggarrange(plot_commodity,plot_forestry,plot_shifting_agr,plot_urban,ncol = 2,nrow = 2)

ggsave('Output//Fig7.jpeg',width=16,height=12,dpi = 400)


#######################################
#####Tables AND FIGURES APPENDIX ######
#######################################

#####Table A.1 ----- 

var_analysis <- df_no_EU %>% filter(year < 2020) %>% 
  dplyr::select(log_deforest_30,pre_EUTR_eu_top_imp_weight,
                log_gdp,log_pop,log_total_export_weight,forest_rents,gdp_growth,
                vdem_polyarchy,vdem_ruleoflaw)

stargazer(var_analysis,type='latex',
          title = 'Description of Variables',
          label = 'Desc_Var',
          covariate.labels = c('Log Deforestation (in ha)',
                               'Dummy EU Trade Dependent',
                               'Log GDP','Log Population','Log Total Timber Exports',
                               'Forest Rents (in \\% of GDP)','GDP Growth (in \\%)',
                               'VDEM Polyarchy Index','VDEM Rule of Law Index'),
          font.size = 'footnotesize', 
          out = 'Output//TableA1.tex')


#####Figure A.1 ----- 
#Select Variables
df_no_EU_pre <- df_no_EU %>% filter(year < 2010)

bal_df <- df_no_EU_pre %>% filter(!is.na(deforest_30)) %>% 
  dplyr::select(pre_EUTR_eu_top_imp_weight, log_gdp ,log_pop ,
                log_total_export_weight , forest_rents , gdp_growth ,
                vdem_polyarchy , vdem_ruleoflaw) %>% drop_na()

#Calculate Sample Means by Treatment Group
pre_match = bal_df  %>% 
  mutate(log_gdp=(log_gdp-mean(log_gdp,na.rm=T))/sd(log_gdp,na.rm = T),
         log_pop=(log_pop-mean(log_pop,na.rm=T))/sd(log_pop,na.rm = T),
         log_total_export_weight=(log_total_export_weight-mean(log_total_export_weight,na.rm=T))/sd(log_total_export_weight,na.rm = T),
         forest_rents=(forest_rents-mean(forest_rents,na.rm=T))/sd(forest_rents,na.rm = T),
         gdp_growth=(gdp_growth-mean(gdp_growth,na.rm=T))/sd(gdp_growth,na.rm = T),
         vdem_polyarchy=(vdem_polyarchy-mean(vdem_polyarchy,na.rm=T))/sd(vdem_polyarchy,na.rm = T),
         vdem_ruleoflaw=(vdem_ruleoflaw-mean(vdem_ruleoflaw,na.rm=T))/sd(vdem_ruleoflaw,na.rm = T))  %>%  
  pivot_longer(
    cols =! "pre_EUTR_eu_top_imp_weight",
    names_to = "Variable",
    values_to = "Value") %>% 
  group_by(pre_EUTR_eu_top_imp_weight, Variable) %>% # Group by the variable name and the treatment status
  dplyr::summarize(Mean=mean(Value), #Step 2: Calculate the mean, Lower and Higher 95% confidence intervals
            LCI=mean(Value)-1.96*plotrix::std.error(Value), 
            HCI=mean(Value)+1.96*plotrix::std.error(Value)) 

#Make Figure
ggplot(data=pre_match) +
  geom_linerange(aes(x=Variable, ymin=LCI, ymax=HCI,color=factor(pre_EUTR_eu_top_imp_weight)), size=1,position = position_dodge(width = 0.1))+
  geom_point(size=2, aes(x=Variable,y=Mean,color=factor(pre_EUTR_eu_top_imp_weight)),position = position_dodge(width = 0.1))+
  coord_flip()+
  scale_color_manual(labels = c("Control Countries","Treated Countries"), 
                     values=c( "#006696","#FFA500")) +
  scale_x_discrete(limits=c('vdem_ruleoflaw','vdem_polyarchy','gdp_growth','forest_rents',
                            "log_total_export_weight","log_pop","log_gdp"),
                   labels=c("log_gdp" = "Log GDP",'log_pop' = 'Log Population' ,"log_total_export_weight" = "Log Total Timber Export",
                            "forest_rents" = "Forest Rents","gdp_growth" = "GDP Growth",
                            "vdem_polyarchy" = "VDM Polyarchy Index","vdem_ruleoflaw" = "VDEM Rule of Law Index")) +
  theme_bw()+
  labs(title="Covariate Balance Pre-EUTR", x="", y="Mean Standarized Value", color="Treatment")

ggsave('Output//FigA1.jpeg',width=7,height=5,dpi=400)



#####Figure A.2 ----- 
avg_share_timber_dest <- df_no_EU %>% mutate(row_share_weight = 100 - eu_share_weight-usa_share_weight-china_share_weight) %>%
  group_by(pre_EUTR_eu_top_imp_weight,year) %>%
  summarise(mean_EU_share = mean(eu_share_weight,na.rm = T),
            mean_USA_share = mean(usa_share_weight,na.rm = T),
            mean_China_share = mean(china_share_weight,na.rm = T),
            mean_ROW_share = mean(row_share_weight,na.rm = T)) %>% drop_na(pre_EUTR_eu_top_imp_weight) %>%
  gather(key = 'Destination', value = 'exp_val',3:6)

ggplot(avg_share_timber_dest,aes(x = year, y = exp_val, color = Destination)) + 
  geom_line(size = 1) + 
  facet_wrap(pre_EUTR_eu_top_imp_weight~.,labeller = as_labeller(c(
    `0` = "Control Countries",
    `1` = "Treated Countries"))
  ) + 
  scale_color_manual(labels = c("China","EU",'Rest of the World','USA'), 
                     values=c( "#D55E00","#0072B2",'#999999','#009E73')) +
  geom_vline(xintercept = 2010, linetype = 'dashed',size = 0.7) + 
  labs(x = 'Year',y = 'Share Timber Exports (in %)') + 
  theme_bw()

ggsave('Output//FigA2.jpeg',width=10,height=6,dpi = 400)

#####Figure B.1 ----- 
ggplot(df_no_EU %>% drop_na(pre_EUTR_eu_top_imp_weight) %>% filter(year < 2010),
       aes(x=year,y=log_deforest_30,color = as.factor(pre_EUTR_eu_top_imp_weight))) + 
  geom_smooth(method = 'loess',se=F) + 
  scale_color_manual(labels = c("Control Countries","Treated Countries"), 
                     values=c( "#006696","#FFA500")) +
  labs(color = 'Treatment Group') + 
  scale_x_continuous(breaks=seq(2001,2009,1)) +
  xlab('Year') + ylab('Log Deforestation') +
  theme_bw()

ggsave('Output//FigB1.jpeg',width=8,height=5,dpi = 400)

#####Figure B.2 ----- 
###Panel A
coef_placebo_ALL <- data.frame(year = seq(2002,2009,1),est = NA, ci.low = NA, ci.high = NA)
treat_seq <- seq(2002,2009,1)

for(i in 1:length(treat_seq)){
  df_placebo <- df_no_EU_pre
  df_placebo$treat_new <- if_else(df_placebo$pre_EUTR_eu_top_imp_weight == 1 & df_placebo$year >= treat_seq[i],1,0)
  reg_res <- felm(log_deforest_30 ~ treat_new |
                    year + iso3c|0|iso3c,df_placebo )
  coef_placebo_ALL[i,2] <- reg_res$coefficients[1]
  coef_placebo_ALL[i,3] <- reg_res$coefficients[1] - 1.96*reg_res$cse[1]
  coef_placebo_ALL[i,4] <- reg_res$coefficients[1] + 1.96*reg_res$cse[1]
}

ggplot(coef_placebo_ALL,aes(x=year,y=est,ymin=ci.low,ymax=ci.high)) + 
  geom_pointrange() + 
  geom_hline(yintercept = 0,linetype = 'dashed') + 
  labs(x = 'Placebo Start Year',y = 'Estimate') + 
  scale_x_continuous(breaks=seq(2002,2009,1)) +
  theme_bw()

ggsave('Output//FigB2a.jpeg',width=7,height=4,dpi = 400)


###Panel B
df_no_EU_pre_2013 <- df_no_EU %>% filter(year < 2013)

coef_placebo_ALL_2013 <- data.frame(year = seq(2002,2013,1),est = NA, ci.low = NA, ci.high = NA)
treat_seq_2013 <- seq(2002,2012,1)

for(i in 1:length(treat_seq_2013)){
  df_placebo <- df_no_EU_pre_2013
  df_placebo$treat_new <- if_else(df_placebo$pre_EUTR_eu_top_imp_weight == 1 & df_placebo$year >= treat_seq_2013[i],1,0)
  reg_res <- felm(log_deforest_30 ~ treat_new |
                    year + iso3c|0|iso3c,df_placebo )
  coef_placebo_ALL_2013[i,2] <- reg_res$coefficients[1]
  coef_placebo_ALL_2013[i,3] <- reg_res$coefficients[1] - 1.96*reg_res$cse[1]
  coef_placebo_ALL_2013[i,4] <- reg_res$coefficients[1] + 1.96*reg_res$cse[1]
}

ggplot(coef_placebo_ALL_2013 %>% filter(year < 2013),aes(x=year,y=est,ymin=ci.low,ymax=ci.high)) + 
  geom_pointrange() + 
  geom_hline(yintercept = 0,linetype = 'dashed') + 
  labs(x = 'Placebo Start Year',y = 'Estimate') + 
  scale_x_continuous(breaks=seq(2002,2012,1)) +
  theme_bw()

ggsave('Output//FigB2b.jpeg',width=7,height=4,dpi = 400)


#####Figure C.1 ----- 
###Panel A
ALL_es_cov <- feols(log_deforest_30  ~ i(year, pre_EUTR_eu_top_imp_weight, 2009) + log_gdp + log_pop + 
                      log_total_export_weight + forest_rents + gdp_growth +
                      vdem_polyarchy + vdem_ruleoflaw|                    ## Other controls
                      year + iso3c,                             ## FEs
                    cluster = ~iso3c,                          ## Clustered SEs
                    data = df_no_EU %>% filter(year < 2020))

jpeg(
  filename = "Output//FigC1a.jpeg",
  width = 10,height = 7,units = "in",res = 400)
iplot(ALL_es_cov,main = '',xlab = 'Year',ylab = 'ATT of EUTR on Log Deforestation',
      pt.join = F,lwd = 1.5,cex =2)
dev.off()

###Panel B
df_no_EU$treat_group <- if_else(df_no_EU$pre_EUTR_eu_top_imp_weight == 1, 2010,0)

ALL_csa_cov <- att_gt(yname = "log_deforest_30",
                      tname = "year",
                      idname = "cown",
                      gname = "treat_group",
                      xformla = ~log_gdp + log_pop + 
                        log_total_export_weight + forest_rents + gdp_growth +
                        vdem_polyarchy + vdem_ruleoflaw,
                      data = df_no_EU %>% filter(year < 2020),
                      est_method = "ipw")

jpeg(
  filename = "Output//FigC1b.jpeg",
  width = 10,height = 7,units = "in",res = 400)
ggdid(ALL_csa_cov,title = "",xlab = 'Year',ylab = 'ATT of EUTR on Log Deforestation',
      xgap = 2)
dev.off()

ALL_csa_cov_agg <- aggte(ALL_csa_cov, type = "simple")
summary(ALL_csa_cov_agg)


###Panel C
df_no_EU <- df_no_EU %>% group_by(iso3c) %>%
  mutate(forest_rents0 = forest_rents[year == 2005][1],
         log_total_export_weight0 = log_total_export_weight[year == 2005][1],
         gdp_growth0 = gdp_growth[year == 2005][1],
         log_gdp0 = log_gdp[year == 2005][1],
         log_pop0 = log_pop[year == 2005][1],
         vdem_polyarchy0 = vdem_polyarchy[year == 2005][1],
         vdem_ruleoflaw0 = vdem_ruleoflaw[year == 2005][1],
         forest_area0 = forest_area[year == 2005][1],
         log_forest_area0 = log_forest_area[year == 2005][1])

ALL_es_cov_pre <- feols(log_deforest_30  ~ i(year, pre_EUTR_eu_top_imp_weight, 2009) + 
                          log_gdp0*year + log_pop0*year + 
                          log_total_export_weight0*year + forest_rents0*year + gdp_growth0*year +
                          vdem_polyarchy0*year + vdem_ruleoflaw0*year |                    ## Other controls
                          year + iso3c,                             ## FEs
                        cluster = ~iso3c,                          ## Clustered SEs
                        data = df_no_EU %>% filter(year < 2020))


jpeg(
  filename = "Output//FigC1c.jpeg",
  width = 10,height = 7,units = "in",res = 400)
iplot(ALL_es_cov_pre,main = '',xlab = 'Year',ylab = 'ATT of EUTR on Log Deforestation',
      pt.join = F,lwd = 1.5,cex =2)
dev.off()

#####Figure C.2 ----- 
###Panel A
ALL_ES_FR_AboveMedian <- feols(log_deforest_30 ~ i(year, pre_EUTR_eu_top_imp_weight, 2009)   |                    ## Other controls
                                 year + iso3c,                             ## FEs
                               cluster = ~iso3c,                          ## Clustered SEs
                               data = df_no_EU %>% filter(year < 2020 & pre_EUTR_forest_rents >= 0.27))

jpeg(
  filename = "Output//FigC2a.jpeg",
  width = 10,height = 7,units = "in",res = 400)
iplot(ALL_ES_FR_AboveMedian,main = '',xlab = 'Year',ylab = 'ATT of EUTR on Log Deforestation',
      pt.join = F,lwd = 1.5,cex =2)
dev.off()

###Panel B
ALL_ES_FR_Top3rd <- feols(log_deforest_30 ~ i(year, pre_EUTR_eu_top_imp_weight, 2009)   |                    ## Other controls
                            year + iso3c,                             ## FEs
                          cluster = ~iso3c,                          ## Clustered SEs
                          data = df_no_EU %>% filter(year < 2020 & pre_EUTR_forest_rents >= 2.22))

jpeg(
  filename = "Output//FigC2b.jpeg",
  width = 10,height = 7,units = "in",res = 400)
iplot(ALL_ES_FR_Top3rd,main = '',xlab = 'Year',ylab = 'ATT of EUTR on Log Deforestation',
      pt.join = F,lwd = 1.5,cex =2)
dev.off()


#####Figure C.3 ----- 
###Panel A
Top3rd_es <- feols(log_deforest_30  ~ i(year, EU_dep_weight_top3rd, 2009)|                    ## Other controls
                     year + iso3c,                             ## FEs
                   cluster = ~iso3c,                          ## Clustered SEs
                   data = df_no_EU %>% filter(year < 2020))

jpeg(
  filename = "Output//FigC3a.jpeg",
  width = 10,height = 7,units = "in",res = 400)
iplot(Top3rd_es,main = 'Top 3rd Percentile EU Trade Share',xlab = 'Year',ylab = 'ATT of EUTR on Log Deforestation',
      pt.join = F,lwd = 1.5,cex =2)
dev.off()

###Panel B
TopQ_es <- feols(log_deforest_30  ~ i(year, EU_dep_weight_topQuart, 2009)|                    ## Other controls
                   year + iso3c,                             ## FEs
                 cluster = ~iso3c,                          ## Clustered SEs
                 data = df_no_EU %>% filter(year < 2020))

jpeg(
  filename = "Output//FigC3b.jpeg",
  width = 10,height = 7,units = "in",res = 400)
iplot(TopQ_es,main = 'Top Quarter Percentile EU Trade Share',xlab = 'Year',ylab = 'ATT of EUTR on Log Deforestation',
      pt.join = F,lwd = 1.5,cex =2)
dev.off()


#####Figure C.4 ----- 
df_no_EU$pre_EUTR_weight_EU_01 <- df_no_EU$pre_EUTR_weight_EU/100

contiuous_es <- feols(log_deforest_30  ~ i(year,pre_EUTR_weight_EU_01 , 2009)|                    ## Other controls
                        year + iso3c,                             ## FEs
                      cluster = ~iso3c,                          ## Clustered SEs
                      data = df_no_EU %>% filter(year < 2020))

jpeg(
  filename = "Output//FigC4.jpeg",
  width = 10,height = 7,units = "in",res = 400)
iplot(contiuous_es,main = 'Continuous EU Trade Share',xlab = 'Year',ylab = 'ATT of EUTR on Log Deforestation',
      pt.join = F,lwd = 1.5,cex =2)
dev.off()

#####Figure C.5 -----
###Panel A
coef_placebo_dist <- data.frame(year = seq(2002,2009,1),est = NA, ci.low = NA, ci.high = NA)
treat_seq <- seq(2002,2009,1)

for(i in 1:length(treat_seq)){
  df_placebo <- df_no_EU_pre
  df_placebo$treat_new <- if_else(df_placebo$year >= treat_seq[i],1,0)
  reg_res <- felm(log_deforest_30 ~ treat_new:log_dist_EU |
                    year + iso3c|0|iso3c,df_placebo )
  coef_placebo_dist[i,2] <- reg_res$coefficients[1]
  coef_placebo_dist[i,3] <- reg_res$coefficients[1] - 1.96*reg_res$cse[1]
  coef_placebo_dist[i,4] <- reg_res$coefficients[1] + 1.96*reg_res$cse[1]
}

ggplot(coef_placebo_dist,aes(x=year,y=est,ymin=ci.low,ymax=ci.high)) + 
  geom_pointrange() + 
  geom_hline(yintercept = 0,linetype = 'dashed') + 
  labs(x = 'Placebo Start Year',y = 'Estimate') + 
  theme_bw()

ggsave('Output//FigC5a.jpeg',width=10,height=6)


###Panel B
geodist_timeFE <- feols(log_deforest_30 ~ i(year, pre_EUTR_eu_top_imp_weight, 2009) + 
                          i(year, log_dist_EU, 2009)|                    ## Other controls
                          year + iso3c,                             ## FEs
                        cluster = ~iso3c,                          ## Clustered SEs
                        data = df_no_EU )

jpeg(
  filename = "Output//FigC5b.jpeg",
  width = 10,height = 7,units = "in",res = 400)
iplot(geodist_timeFE,main = '',xlab = 'Year',ylab = 'ATT of EUTR on Log Deforestation',
      pt.join = F,lwd = 1.5,cex =2)
dev.off()

###Panel C
geodist_timetrend <- feols(log_deforest_30 ~ i(year, pre_EUTR_eu_top_imp_weight, 2009) + 
                             year*log_dist_EU|                    ## Other controls
                             year + iso3c,                             ## FEs
                           cluster = ~iso3c,                          ## Clustered SEs
                           data = df_no_EU )

jpeg(
  filename = "Output//FigC5c.jpeg",
  width = 10,height = 7,units = "in",res = 400)
iplot(geodist_timetrend,main = '',xlab = 'Year',ylab = 'ATT of EUTR on Log Deforestation',
      pt.join = F,lwd = 1.5,cex =2)
dev.off()

#####Figure C.6 ----- 
###Panel A
ALL_csa <- att_gt(yname = "log_deforest_30",
                  tname = "year",
                  idname = "cown",
                  gname = "treat_group",
                  base_period = 'universal',
                  xformla = ~1,
                  data = df_no_EU %>% filter(year < 2020),
                  est_method = "reg")

jpeg(
  filename = "Output//FigC6a.jpeg",
  width = 10,height = 7,units = "in",res = 400)
ggdid(ALL_csa,title = "",xlab = 'Year',ylab = 'ATT of EUTR on Log Deforestation',
      xgap = 2)
dev.off()

###Panel B
defor_fect <- fect(log_deforest_30 ~ treat_approved_weight_topImp,
                   data = df_no_EU %>% filter(year < 2020) %>% distinct(iso3c,year,.keep_all = T), index = c("iso3c","year"), na.rm = T,
                   method = "fe", force = "two-way", se = TRUE, parallel = TRUE, nboots = 200)

jpeg(
  filename = "Output//FigC6b.jpeg",
  width = 10,height = 7,units = "in",res = 400)
plot(defor_fect, main = "", ylab = "ATT of EUTR on Log Deforestation", 
     cex.main = 0.8, cex.lab = 0.8, cex.axis = 0.8,stats = "F.p")
dev.off()

#####Figure C.7 ----- 
na_ids <- unique(df_no_EU$iso3c[is.na(df_no_EU$log_deforest_30) | is.na(df_no_EU$treat_approved_weight_topImp)])
count_obs <- df_no_EU %>% filter(!(iso3c %in% na_ids)) %>% filter(year < 2020) %>% group_by(iso3c) %>% count() %>% ungroup() %>%
  filter(n < 19) %>% dplyr::select(iso3c) %>% as.matrix() %>% as.vector()
df_synth_did <- df_no_EU %>% filter(!(iso3c %in% c(na_ids,count_obs))) %>% filter(year < 2020) %>% 
  distinct(year,iso3c,.keep_all = T) %>% as.data.frame()

setup <-  panel.matrices(df_synth_did,unit = 'iso3c',time = 'year',outcome = 'log_deforest_30',treatment = 'treat_approved_weight_topImp')
tau.hat = synthdid_estimate(setup$Y, setup$N0, setup$T0)
se = sqrt(vcov(tau.hat, method='jackknife'))
sprintf('point estimate: %1.2f', tau.hat)
sprintf('95%% CI (%1.2f, %1.2f)', tau.hat - 1.96 * se, tau.hat + 1.96 * se)

jpeg(
  filename = "Output//FigC7.jpeg",
  width = 10,height = 7,units = "in",res = 400)
synthdid_plot(tau.hat,trajectory.linetype = 1, line.width=1.25,
              control.name='Control', treated.name='Treated')  + theme_bw()
dev.off()

#####Figure C.8 ----- 
###Panel A
defor50_es <- feols(log_deforest_50  ~ i(year, pre_EUTR_eu_top_imp_weight, 2009) |                    ## Other controls
                      year + iso3c,                             ## FEs
                    cluster = ~iso3c,                          ## Clustered SEs
                    data = df_no_EU)

jpeg(
  filename = "Output//FigC8a.jpeg",
  width = 10,height = 7,units = "in",res = 400)
iplot(defor50_es,main = '',xlab = 'Year',ylab = 'ATT of EUTR on Log Deforestation',
      pt.join = F,lwd = 1.5,cex =2)
dev.off()

###Panel B
defor75_es <- feols(log_deforest_75  ~ i(year, pre_EUTR_eu_top_imp_weight, 2009) |                    ## Other controls
                      year + iso3c,                             ## FEs
                    cluster = ~iso3c,                          ## Clustered SEs
                    data = df_no_EU)

jpeg(
  filename = "Output//FigC8b.jpeg",
  width = 10,height = 7,units = "in",res = 400)
iplot(defor75_es,main = '',xlab = 'Year',ylab = 'ATT of EUTR on Log Deforestation',
      pt.join = F,lwd = 1.5,cex =2)
dev.off()

#####Figure C.9 ----- 
###Panel A
area_es_1 <- feols(log_deforest_percent_30 ~ i(year, pre_EUTR_eu_top_imp_weight, 2009)|                    ## Other controls
                     year + iso3c,      
                   data = df_no_EU %>% filter(year < 2020),
                   cluster = "iso3c") 

jpeg(
  filename = "Output//FigC9a.jpeg",
  width = 10,height = 7,units = "in",res = 400)
iplot(area_es_1,main = '',xlab = 'Year',ylab = 'ATT of EUTR on Percentage of Deforested Area',
      pt.join = F,lwd = 1.5,cex =2)
dev.off()

###Panel B
area_es_2 <- feols(log_deforest_percent_50 ~ i(year, pre_EUTR_eu_top_imp_weight, 2009)|                    ## Other controls
                     year + iso3c,      
                   data = df_no_EU %>% filter(year < 2020),
                   cluster = "iso3c") 

jpeg(
  filename = "Output//FigC9b.jpeg",
  width = 10,height = 7,units = "in",res = 400)
iplot(area_es_2,main = '',xlab = 'Year',ylab = 'ATT of EUTR on Percentage of Deforested Area',
      pt.join = F,lwd = 1.5,cex =2)
dev.off()

#####Figure C.10 -----
prime_deforest_model <- feols(log(prime_deforest+1) ~ i(year, pre_EUTR_eu_top_imp_weight, 2009) |                    ## Other controls
                                year + iso3c,                             ## FEs
                              cluster = ~iso3c,                          ## Clustered SEs
                              data = df_no_EU %>% filter(year < 2020))

jpeg(
  filename = "Output//FigC10.jpeg",
  width = 10,height = 7,units = "in",res = 400)
iplot(prime_deforest_model,main = 'The Effect of the EUTR on Primary Forest Loss',xlab = 'Year',ylab = 'ATT of EUTR on Log Deforestation',
      pt.join = F,lwd = 1.5,cex =2)
dev.off()

#####Figure C.11 -----
###Panel A
no_vpaall_model <- feols(log_deforest_30  ~ i(year, pre_EUTR_eu_top_imp_weight, 2009) |                    ## Other controls
                           year + iso3c,                             ## FEs
                         cluster = ~iso3c,                          ## Clustered SEs
                         data = df_no_EU %>% filter(year < 2020) %>%
                           filter(vpa_country != 1))

jpeg(
  filename = "Output//FigC11a.jpeg",
  width = 10,height = 7,units = "in",res = 400)
iplot(no_vpaall_model,main = 'The Effect of the EUTR on Deforestation (No VPA Country)',xlab = 'Year',ylab = 'ATT of EUTR on Log Deforestation',
      pt.join = F,lwd = 1.5,cex =2)
dev.off()

#Panel B
eea_countries <- c("ISL",'LIE',"NOR")

no_eea_model <- feols(log_deforest_30  ~ i(year, pre_EUTR_eu_top_imp_weight, 2009) |                    ## Other controls
                        year + iso3c,                             ## FEs
                      cluster = ~iso3c,                          ## Clustered SEs
                      data = df_no_EU %>% filter(year < 2020) %>%
                        filter(!(iso3c %in% eea_countries)))

jpeg(
  filename = "Output//FigC11b.jpeg",
  width = 10,height = 7,units = "in",res = 400)
iplot(no_eea_model,main = 'The Effect of the EUTR on Deforestation (No EEA Countries)',xlab = 'Year',ylab = 'ATT of EUTR on Log Deforestation',
      pt.join = F,lwd = 1.5,cex =2)
dev.off()

###Panel C
oecd_non_EU <- c('AUS','CAN','CHL','COL','CRI','ISL','ISR','JPN','KOR','MEX','NZL','NOR','CHE','TUR','USA')

no_oecd_model <- feols(log_deforest_30  ~ i(year, pre_EUTR_eu_top_imp_weight, 2009) |                    ## Other controls
                         year + iso3c,                             ## FEs
                       cluster = ~iso3c,                          ## Clustered SEs
                       data = df_no_EU %>% filter(year < 2020) %>%
                         filter(!(iso3c %in% oecd_non_EU)))

jpeg(
  filename = "Output//FigC11c.jpeg",
  width = 10,height = 7,units = "in",res = 400)
iplot(no_oecd_model,main = 'The Effect of the EUTR on Deforestation (No OECD Countries)',xlab = 'Year',ylab = 'ATT of EUTR on Log Deforestation',
      pt.join = F,lwd = 1.5,cex =2)
dev.off()


#####Figure C.12 -----
ALL_es_km2 <- feols(log_deforest_30_km2  ~ i(year, pre_EUTR_eu_top_imp_weight, 2009) |                    ## Other controls
                      year + iso3c,                             ## FEs
                    cluster = ~iso3c,                          ## Clustered SEs
                    data = df_no_EU)

jpeg(
  filename = "Output//FigC12.jpeg",
  width = 10,height = 7,units = "in",res = 400)
iplot(ALL_es_km2,main = '',xlab = 'Year',ylab = 'ATT of EUTR on Log Deforestation',
      pt.join = F,lwd = 1.5,cex =2)
dev.off()

#####Figure C.13 -----
###Panel A

ALL_es_cov_pta <- feols(log_deforest_30  ~ i(year, pre_EUTR_eu_top_imp_weight, 2009) + log_gdp + log_pop + 
                          log_total_export_weight + forest_rents + gdp_growth +
                          vdem_polyarchy + vdem_ruleoflaw + num_ptas_signed|                    ## Other controls
                          year + iso3c,                             ## FEs
                        cluster = ~iso3c,                          ## Clustered SEs
                        data = df_no_EU %>% filter(year < 2020))

jpeg(
  filename = "Output//FigC13a.jpeg",
  width = 10,height = 7,units = "in",res = 400)
iplot(ALL_es_cov_pta,main = '',xlab = 'Year',ylab = 'ATT of EUTR on Log Deforestation',
      pt.join = F,lwd = 1.5,cex =2)
dev.off()

###Panel B
ALL_es_cov_pta_env <- feols(log_deforest_30  ~ i(year, pre_EUTR_eu_top_imp_weight, 2009) + log_gdp + log_pop + 
                              log_total_export_weight + forest_rents + gdp_growth +
                              vdem_polyarchy + vdem_ruleoflaw + num_ptas_signed_env|                    ## Other controls
                              year + iso3c,                             ## FEs
                            cluster = ~iso3c,                          ## Clustered SEs
                            data = df_no_EU %>% filter(year < 2017))

jpeg(
  filename = "Output//FigC13b.jpeg",
  width = 10,height = 7,units = "in",res = 400)
iplot(ALL_es_cov_pta_env,main = '',xlab = 'Year',ylab = 'ATT of EUTR on Log Deforestation',
      pt.join = F,lwd = 1.5,cex =2)
dev.off()

###Panel C

ALL_es_cov_pta_envrisk <- feols(log_deforest_30  ~ i(year, pre_EUTR_eu_top_imp_weight, 2009) + log_gdp + log_pop + 
                                  log_total_export_weight + forest_rents + gdp_growth +
                                  vdem_polyarchy + vdem_ruleoflaw + num_ptas_signed_envrisk|                    ## Other controls
                                  year + iso3c,                             ## FEs
                                cluster = ~iso3c,                          ## Clustered SEs
                                data = df_no_EU %>% filter(year < 2017))

jpeg(
  filename = "Output//FigC13c.jpeg",
  width = 10,height = 7,units = "in",res = 400)
iplot(ALL_es_cov_pta_envrisk,main = '',xlab = 'Year',ylab = 'ATT of EUTR on Log Deforestation',
      pt.join = F,lwd = 1.5,cex =2)
dev.off()

###Panel D

ALL_es_cov_pta_envforest <- feols(log_deforest_30  ~ i(year, pre_EUTR_eu_top_imp_weight, 2009) + log_gdp + log_pop + 
                                    log_total_export_weight + forest_rents + gdp_growth +
                                    vdem_polyarchy + vdem_ruleoflaw + num_ptas_signed_envforest|                    ## Other controls
                                    year + iso3c,                             ## FEs
                                  cluster = ~iso3c,                          ## Clustered SEs
                                  data = df_no_EU %>% filter(year < 2017))

jpeg(
  filename = "Output//FigC13d.jpeg",
  width = 10,height = 7,units = "in",res = 400)
iplot(ALL_es_cov_pta_envforest,main = '',xlab = 'Year',ylab = 'ATT of EUTR on Log Deforestation',
      pt.join = F,lwd = 1.5,cex =2)
dev.off()

#####Figure D.1 -----
dyad_df <- read_csv('dyad_timber_trade_data.csv')

eu_imp_year <- dyad_df %>% filter(imp_cown == 999) %>% group_by(year) %>%
  summarise(weight = sum(total_export_weight,na.rm = T),
            value = sum(total_export_value,na.rm = T)) %>%
  mutate(weight_scale = weight/116281555*100) %>%
  mutate(name = 'EU Timber Imports')

total_imp_year <- dyad_df %>% filter(imp_cown != 999) %>% group_by(year) %>%
  summarise(weight = sum(total_export_weight,na.rm = T),
            value = sum(total_export_value,na.rm = T)) %>%
  mutate(weight_scale = weight/251073317*100) %>%
  mutate(name = 'Non-EU Timber Imports')

data_trade_dev <- rbind(total_imp_year,eu_imp_year)

ggplot(data_trade_dev %>% filter(year < 2020),
       aes(x = year,y = weight_scale,shape = name,linetype = name)) + 
  geom_line(size =1) + 
  geom_point(size = 3) + 
  labs(linetype = '' , shape = '',x = 'Year', y = 'Standardized Import Growth (Reference Year 2009)') +
  theme_bw()

ggsave('Output//FigD1.jpeg',width=12,height=6,dpi = 400)


#####Figure D.2 -----
###Panel A
wood_imports_EU <- read_excel('wood_imports_EU.xls') 

wi_EU_df <- wood_imports_EU %>% 
  mutate(Unit = if_else(Unit == '1000 USD','value','weight')) %>%
  pivot_wider(names_from = Unit,values_from = Value,id_cols = c(Area,Year,Item)) %>%
  mutate(price = value/weight,
         Year = as.numeric(Year))  %>%
  group_by(Year,Item)%>%
  mutate(total_value = sum(value,na.rm = T),
         share = value/total_value) %>%
  dplyr::summarize(avg_price = mean(price*share,na.rm = T)) %>% ungroup() %>%
  mutate(avg_price_scale = if_else(Item == 'Roundwood', avg_price/0.003195949*100,avg_price/0.011442675*100)) 

ggplot(wi_EU_df, aes(Year, avg_price_scale, colour = Item)) +
  geom_line(linewidth = 1) + 
  geom_vline(xintercept = 2010,size = 0.75,linetype = 'dashed') + 
  scale_x_continuous(breaks = seq(min(wi_EU_df$Year), max(wi_EU_df$Year), 1)) +
  scale_colour_manual(values = c("#0072B2", "#D55E00")) +
  labs(colour = 'Timber Type',x='Year',y = 'Weighted Average Price Index (2010 = 100)') +
  theme_bw(base_size = 12) +
  theme(panel.grid.minor = element_blank(),
        axis.text.x = element_text(angle = 45, hjust = 1),
        legend.position = "bottom")

ggsave('Output//FigD1a.jpeg',width=10,height=6,dpi = 400)

###Panel B
eu_producer_cpi <- tibble(cpi = c(89.1,92.2,94.7,94.7,94.4	,95.4,	96.0	,95.7,	95.3,	95.8,	95.4,	95.0,	94.8,	96.2,	97.2,	98.6,	99.5	,100.4	,100.6,	100.3,	100.3,	100.4,	100.0,	99.3	,98.4,	98.6,	99.1,	99.5,	99.7,	100.5	,101.1	,102.1,	103.6	,106.0	,107.6	,107.8	,107.4	,106.4,	104.4	,103.0),
                          q_seq = seq(from = as.yearqtr("2010 Q1"),
                                      to   = as.yearqtr("2019 Q4"),
                                      by   = 1/4))

ggplot(eu_producer_cpi,aes(x = q_seq,y=cpi)) + 
  geom_line(size = 1) + 
  scale_x_continuous(
    breaks = as.numeric(eu_producer_cpi$q_seq),
    labels = function(x) format(as.yearqtr(x), "%Y Q%q"),
    expand = expansion(mult = c(0.01, 0.01))
  ) +
  labs(x = 'Quarters',y = 'Producer Price Index (2015 = 1)') +
  theme_bw() +
  theme(
    axis.text.x = element_text(angle = 45, hjust = 1),
    panel.grid.minor = element_blank(),
    legend.position = "bottom"
  )

ggsave('Output//FigD1b.jpeg',width=10,height=6,dpi = 400)


###Panel C
imp_eu_df <- dyad_df %>% filter(in_EU_imp == 1 & in_EU_exp ==0) %>% group_by(year) %>%
  summarise(total_weight = sum(total_export_weight),
            total_value = sum(total_export_value)) %>% 
  mutate(price = total_value / total_weight) %>% ungroup() %>%
  group_by(year) %>%
  mutate(weight_scale = total_weight/41768840*100,
         price_scale = price/0.4171156*100,
         Type = 'Imports to EU') 

intra_eu_df <- dyad_df %>% filter(in_EU_exp ==1 & in_EU_imp == 1) %>% group_by(year) %>%
  summarise(total_weight = sum(total_export_weight),
            total_value = sum(total_export_value)) %>% 
  mutate(price = total_value / total_weight) %>%
  group_by(year) %>%
  mutate(weight_scale = total_weight/95467770*100,
         price_scale = price/0.3544334*100,
         Type = 'Intra-EU Imports') 

ggplot(rbind(imp_eu_df,intra_eu_df),aes(x = year,y = weight_scale,colour = Type))  +
  geom_line(linewidth = 1) + 
  geom_vline(xintercept = 2010,size = 0.75,linetype = 'dashed') + 
  scale_x_continuous(breaks = seq(2001,2020, 1)) +
  scale_colour_manual(values = c("#0072B2", "#D55E00")) +
  labs(colour = 'Trade Type',x='Year',y = 'Timber Import Index (2010 = 100)') +
  theme_bw(base_size = 12) +
  theme(panel.grid.minor = element_blank(),
        axis.text.x = element_text(angle = 45, hjust = 1),
        legend.position = "bottom")

ggsave('Output//FigD1c.jpeg',width=10,height=6,dpi = 400)

#####Figure D.3-----
##Total Exports
total_exports <- feols(log_total_export_weight  ~ i(year, pre_EUTR_eu_top_imp_weight, 2009) |                    ## Other controls
                         year + iso3c,                             ## FEs
                       cluster = ~iso3c,                          ## Clustered SEs
                       data = df_no_EU %>% filter(year < 2020))

agg_totalexports <- aggr_es(total_exports,period = 'post',aggregation = 'mean')


##GDP Growth 
gdp_growth <- feols(gdp_growth  ~ i(year, pre_EUTR_eu_top_imp_weight, 2009) |                    ## Other controls
                      year + iso3c,                             ## FEs
                    cluster = ~iso3c,                          ## Clustered SEs
                    data = df_no_EU %>% filter(year < 2020))

agg_gdp_growth <- aggr_es(gdp_growth,period = 'post',aggregation = 'mean')

##Timber Production 
forest_production_cubic <- feols(log_forestery_production_cubic  ~ i(year, pre_EUTR_eu_top_imp_weight, 2009) |                    ## Other controls
                                   year + iso3c,                             ## FEs
                                 cluster = ~iso3c,                          ## Clustered SEs
                                 data = df_no_EU %>% filter(year < 2020))

agg_forest_production_cubic <- aggr_es(forest_production_cubic,period = 'post',aggregation = 'mean')

forest_production_tons <- feols(log_forestery_production_tons  ~ i(year, pre_EUTR_eu_top_imp_weight, 2009) |                    ## Other controls
                                  year + iso3c,                             ## FEs
                                cluster = ~iso3c,                          ## Clustered SEs
                                data = df_no_EU %>% filter(year < 2020))

agg_forest_production_tons <- aggr_es(forest_production_tons,period = 'post',aggregation = 'mean')


##Wildfires 
forest_fire <- feols(log(burned_forest_area+1)  ~ i(year, pre_EUTR_eu_top_imp_weight, 2009) |                    ## Other controls
                       year + iso3c,                             ## FEs
                     cluster = ~iso3c,                          ## Clustered SEs
                     data = df_no_EU %>% filter(year < 2020))

agg_forest_fire <- aggr_es(forest_fire,period = 'post',aggregation = 'mean')


#### Make Plot 
alt_exp_df <- dplyr::bind_rows(agg_gdp_growth,agg_forest_production_cubic,agg_forest_production_tons,agg_totalexports,agg_forest_fire) %>%
  mutate(DV = c('GDP Growth','Log Forest Goods Production (in m3)','Log Forest Goods Production (in t)',
                'Log Total Timber Exports (in t)', 'Log Forest Area Burned (in ha)'))

ggplot(alt_exp_df,aes(x = DV, ymin = conf.low,ymax = conf.high,y=estimate)) + 
  geom_pointrange(size = 0.8) + 
  geom_hline(yintercept = 0,linetype = 'dashed') + 
  ylab('Effect Size') + xlab('') +
  theme_bw() + coord_flip()

ggsave('Output//FigD3.jpeg',width=8,height=4,dpi = 400)

#####Figure E.1 -----
df_no_EU_2019 <- df_no_EU_2019 %>%
  group_by(iso3c) %>%   # replace with your actual unit identifier
  mutate(china_share_weight0 = ifelse(
    any(year == 2008),
    china_share_weight[year == 2008],
    NA_real_
  )) %>%
  ungroup()


###Panel A
china_share_change <- feols(china_share_weight  ~ i(year, pre_EUTR_eu_top_imp_weight, 2009)  |                    ## Other controls
                              year + iso3c,                             ## FEs
                            cluster = ~iso3c,                          ## Clustered SEs
                            data = df_no_EU_2019)


jpeg(
  filename = "Output//FigE1a.jpeg",
  width = 10,height = 7,units = "in",res = 400)
iplot(china_share_change,main = 'Effect of EUTR on Timber Export Share to China',xlab = 'Year',ylab = 'Difference Timber Export Share to China',
      pt.join = F,lwd = 1.5,cex =2)
dev.off()

###Panel B
china_share_change_rob <- feols(china_share_weight  ~ i(year, pre_EUTR_eu_top_imp_weight, 2009) + 
                                  i(year, china_share_weight0, 2009) |                    ## Other controls
                                  year + iso3c,                             ## FEs
                                cluster = ~iso3c,                          ## Clustered SEs
                                data = df_no_EU_2019)

jpeg(
  filename = "Output//FigE1b.jpeg",
  width = 10,height = 7,units = "in",res = 400)
iplot(china_share_change_rob,main = 'Effect of EUTR on Timber Export Share to China',xlab = 'Year',ylab = 'Difference Timber Export Share to China',
      pt.join = F,lwd = 1.5,cex =2)
dev.off()

#####Figure E.2 -----
###Panel A
#Run Model
china_share_change_het1_rob <- feols(china_share_weight ~ i(year, pre_EUTR_eu_top_imp_weight, 2009)*pre_EUTR_vdem_pol_corrupt + 
                                       i(year, china_share_weight0, 2009)|                    ## Other controls
                                       year + iso3c,                             ## FEs
                                     cluster = ~iso3c,                          ## Clustered SEs
                                     data = df_no_EU_2019)

#Make Plot
het1_plot_chn_rob <- plot_slopes(
  china_share_change_het1_rob,
  variables = 'pre_EUTR_eu_top_imp_weight',
  condition = list(
    year = unique(df_no_EU_2019$year),
    pre_EUTR_vdem_pol_corrupt = c(0.86, 0.65, 0.2)  # Now: High → Medium → Low
  )
) +
  geom_hline(yintercept = 0, linetype = 'dashed') +
  labs(
    title = 'A) Political Corruption',
    x = 'Year',
    y = 'Difference in Share of Timber Exports to China'
  ) +
  theme_bw() +
  ylim(c(-25,40)) +
  theme(
    legend.position = "none",
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank()
  ) +
  scale_color_manual(values = corruption_colors) +
  scale_x_discrete(breaks = seq(min(df_no_EU_2019$year), max(df_no_EU_2019$year), by = 3)) +
  facet_wrap(
    ~ factor(pre_EUTR_vdem_pol_corrupt, levels = c(0.86, 0.65, 0.2)),  # Force order inside plotting
    ncol = 3,
    labeller = as_labeller(c(
      `0.86` = "High corruption score (0.86)",
      `0.65` = "Median corruption score (0.65)",
      `0.2`  = "Low corruption score (0.2)"
    ))
  )

#Pooled Results
avg_slopes(china_share_change_het1_rob, 
                               variables = "pre_EUTR_eu_top_imp_weight", 
                               newdata = datagrid(
                                 year = seq(2010,2019,1),
                                 pre_EUTR_vdem_pol_corrupt = c(0.2, 0.65, 0.86)
                               ),by = "pre_EUTR_vdem_pol_corrupt")

###Panel B
#Run Model
china_share_change_het2_rob <- feols(china_share_weight ~ i(year, pre_EUTR_eu_top_imp_weight, 2009)*pre_EUTR_vdem_ruleoflaw  + 
                                       i(year, china_share_weight0, 2009) |                    ## Other controls
                                       year + iso3c,                             ## FEs
                                     cluster = ~iso3c,                          ## Clustered SEs
                                     data = df_no_EU_2019)

#Make Figure
het2_plot_chn_rob <- plot_slopes(china_share_change_het2_rob,variables = 'pre_EUTR_eu_top_imp_weight' , condition = list(
  year = unique(df_no_EU_2019$year),  # All unique years from the dataset
  pre_EUTR_vdem_ruleoflaw = c(0.13, 0.44, 0.83)  # Values of interest
)) + geom_hline(yintercept = 0, linetype = 'dashed') +
  labs(title = 'B) Rule of Law' ,x = 'Year', y = 'Difference in Share of Timber Exports to China') +
  theme_bw()  +ylim(c(-25,40)) +
  theme(legend.position = "none") +  # Move legend for better readability
  scale_color_manual(values = rol_colors) +  # Adjust colors for clarity
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) + 
  scale_x_discrete(breaks = seq(min(df_no_EU_2019$year), max(df_no_EU_2019$year), by = 3)) +  # Show only every third year
  facet_wrap(~pre_EUTR_vdem_ruleoflaw, ncol = 3, labeller = as_labeller(c(
    `0.13` = "Low Rule of Law score (0.13)",
    `0.44` = "Median Rule of Law score (0.44)",
    `0.83` = "High Rule of Law score (0.83)"
  )))  # Separates estimates into clear facets

#Pooled Results
avg_slopes(china_share_change_het2_rob, 
                                 variables = "pre_EUTR_eu_top_imp_weight", 
                                 newdata = datagrid(
                                   year = seq(2010,2019,1),
                                   pre_EUTR_vdem_ruleoflaw = c(0.13, 0.44, 0.83)
                                 ),by = "pre_EUTR_vdem_ruleoflaw")

###Panel C
#Run Model
china_share_change_het3_rob <- feols(china_share_weight ~ i(year, pre_EUTR_eu_top_imp_weight, 2009)*pre_EUTR_state_cap   + 
                                       i(year, china_share_weight0, 2009)|                    ## Other controls
                                       year + iso3c,                             ## FEs
                                     cluster = ~iso3c,                          ## Clustered SEs
                                     data = df_no_EU_2019)

#Make Figure
het3_plot_china_rob <- plot_slopes(china_share_change_het3_rob,variables = 'pre_EUTR_eu_top_imp_weight' , condition = list(
  year = unique(df_no_EU_2019$year),  # All unique years from the dataset
  pre_EUTR_state_cap = c(-0.44, 0.23, 0.9)  # Values of interest
)) + geom_hline(yintercept = 0, linetype = 'dashed') +
  labs(title = 'C) State Capacity' ,x = 'Year', y = 'Difference in Share of Timber Exports to China') +
  theme_bw()  + ylim(c(-25,40)) +
  theme(legend.position = "none") +  # Move legend for better readability
  scale_color_manual(values = rol_colors) +  # Adjust colors for clarity
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) + 
  scale_x_discrete(breaks = seq(min(df_no_EU_2019$year), max(df_no_EU_2019$year), by = 3)) +  # Show only every third year
  facet_wrap(~pre_EUTR_state_cap, ncol = 3, labeller = as_labeller(c(
    `-0.44` = "Low State Capacity score (-0.44)",
    `0.23` = "Median State Capacity score (0.23)",
    `0.9` = "High State Capacity score (0.9)"
  )))  # Separates estimates into clear facets

#Pooled Result
avg_slopes(china_share_change_het3_rob, 
                                 variables = "pre_EUTR_eu_top_imp_weight", 
                                 newdata = datagrid(
                                   year = seq(2010,2019,1),
                                   pre_EUTR_state_cap = c(-0.44, 0.23, 0.9)
                                 ),by = "pre_EUTR_state_cap")

#Combine Figures
ggpubr::ggarrange(het1_plot_chn_rob,het2_plot_chn_rob,het3_plot_china_rob,nrow = 3,ncol = 1)

ggsave('Output//FigE2.jpeg',width = 18,height = 18,dpi = 400)


